# -*- coding: utf-8 -*-
"""
Created on Thu Jun 12 11:44:42 2014

@author: Maxim
"""

from perf_tools import FlightMechanics, ISAtmosphere, FlightConditions
from scipy.optimize import fsolve, minimize_scalar, brentq
import numpy as np
import constants

class PerformanceResults:
    def __init__(self):
        self.velocity              = 0.0
        self.Mach                  = 0.0
        self.SAR                   = 0.0
        self.TAS                   = 0.0
        self.EAS                   = 0.0
        self.altitude              = 0.0
        self.density               = 0.0
        self.climbRate             = 0.0
        self.climbAngle            = 0.0
        self.fuelFlow              = 0.0
        self.TSFC                  = 0.0
        self.thrust                = 0.0
        self.lift                  = 0.0
        self.drag                  = 0.0
        self.LD                    = 0.0
        self.turnRadius            = 0.0
        self.bankAngle             = 0.0
        self.loadFactor            = 1.0
        self.turnRate              = 0.0
        self.weight                = 0.0
        self.emptyWeight           = 0.0
        self.fuelWeight            = 0.0
        self.fuelWeightFraction    = 0.0
        self.payloadWeight         = 0.0
        self.reynoldsNumber        = 0.0
        self.liftCoefficient       = 0.0
        self.parasiteDrag          = 0.0
        self.induceedDrag          = 0.0
        self.dragCoefficient       = 0.0
        self.parasiteDragCoeffient = 0.0
        self.induceedDragCoeefient = 0.0


    def process_data(self):
        self.EAS = self.velocity * (self.density/1.2255)**0.5

    def __repr__1(self):
        out = 'performance data\n================\n'
        for attr, value in self.__dict__.iteritems():
            if type(value) is float or type(value) is np.float64:
                out += '{:<15} = {:<12.5f}\n'.format(attr,value)
        return out

    def __repr__(self):
        self.process_data()
        out = 'performance data\n================\n'
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('altitude',self.altitude,'m')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('density',self.density,'kg/m3')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('TAS',self.velocity,'m/s')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('EAS',self.EAS,'m/s')
        #out += '{:<23} = {:<12.2f}{:<6}\n'.format('IAS',self.IAS,'m/s')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('Mach',self.Mach,'')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('L/D',self.LD,'')
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('lift',self.lift,'N')
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('drag',self.drag,'N')
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('thrust',self.thrust,'N')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('SAR',self.SAR,'m/kg')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('fuel flow',self.fuelFlow,'kg/s')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('TSFC',self.TSFC,'...')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('climb rate',self.climbRate,'m/s')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('climb angle',self.climbAngle,'deg')
        #out += '{:<23} = {:<12.5f}{:>6}\n'.format('turn radius',self.SAR,'m/kg')
        #out += '{:<23} = {:<12.5f}{:>6}\n'.format('bank angle',self.SAR,'m/kg')
        #out += '{:<23} = {:<12.5f}{:>6}\n'.format('turn rate',self.SAR,'m/kg')
        #out += '{:<23} = {:<12.5f}{:>6}\n'.format('load factor',self.loadFactor,'')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('weight',self.weight,'kg')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('empty weight',self.emptyWeight,'kg')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('fuel weight',self.fuelWeight,'kg')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('fuel weight fraction',self.fuelWeightFraction,'')
        out += '{:<23} = {:<12.2f}{:<6}\n'.format('payload weight',self.payloadWeight,'kg')
        out += '{:<23} = {:<12.2e}{:<6}\n'.format('reynolds number',self.reynoldsNumber,'')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('lift coefficient',self.liftCoefficient,'')
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('parasite drag',self.parasiteDrag,'N')  
        out += '{:<23} = {:<12.1f}{:<6}\n'.format('induceed drag',self.induceedDrag,'N')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('drag coefficient',self.dragCoefficient,'')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('parasite drag coeffient',self.parasiteDragCoeffient,'')
        out += '{:<23} = {:<12.4f}{:<6}\n'.format('induceed drag coeefient',self.induceedDragCoeefient,'')
        return out
    def display(self):
        print self.__repr__() 



class SteadyLevelFlight(FlightMechanics):
    """
    Set of methods for aircraft performance analysis in steady level flight. 
    
    Parameters
    ----------
    
    aircraft : object
        pre-loaded aircraft object.
    CLmax : float
        maximum 3D lift coefficient. This value is used for stall speed calculation.
    
    Examples
    --------
    
    Load aircraft
    
    >>> ac = aircraft_FW.load('X45C')
    
    Initialize steady level flight class

    >>> CLmax3D = 1.6
    >>> slf = SteadyLevelFlight(ac, CLmax3D)
    >>> altitude = 10000 # 10km
    """
    def run_max_TAS(self,altitude):
        """
        Calculated maximum true airspeed at given altitude: speed at which 
        maximum thrust at given altitude is equal to drag. 
        
        Parameters
        ----------
        
        altitude : float, m
        
        Examples
        --------
        
        Get conditions at which true airspeed is maximized.
        
        >>> slf.run_max_TAS(altitude)
        performance data
        =========
        LD              = 4.56022     
        turnRadius      = 0.00000     
        TSFC            = 0.00033     
        density         = 0.41270     
        fuelFlow        = 0.55606     
        loadFactor      = 1.00000     
        drag            = 16466.27101 
        lift            = 75089.77219 
        SAR             = 495.72897   
        thrust          = 16466.27101 
        velocity        = 275.65341   
        EAS             = 0.00000     
        turnRate        = 0.00000     
        climbRate       = 0.00000     
        climbAngle      = 0.00000     
        bankAngle       = 0.00000     
        Mach            = 0.92051     
        
        """
        atm = ISAtmosphere(altitude)
        def velocity_eqn(V,altitude):
            Treq = self.get_required_thrust(V,altitude)
            Tav = self.tm.get_thrust_available(altitude)
            return Treq-Tav
        V0 = atm.soundSpeed*0.3
        V2 = atm.soundSpeed
        try:
            V = brentq(velocity_eqn,V0,V2,xtol=1e-3,args=(altitude,))
        except ValueError:
            V = fsolve(velocity_eqn,atm.soundSpeed*0.9,(altitude,))[0]
        self.set_results(V,altitude)
        self.results.thrust = self.get_required_thrust(V,altitude)
        return self.results
    
    def run_min_TAS(self,altitude):
        atm = ISAtmosphere(altitude)
        def velocity_eqn(V,altitude):
            Treq = self.get_required_thrust(V,altitude)
            Tav = self.tm.get_thrust_available(altitude)
            return Treq-Tav
        V0 = atm.soundSpeed*0.3
        V = fsolve(velocity_eqn,V0,(altitude,))[0]
        self.set_results(V,altitude)
        self.results.thrust = self.tm.totalThrust
        return self.results
    
    def run_max_SAR(self,altitude):
        """
        Calculate flight conditions at which specific air range is maximized:
        maximum range condition
        
        Parameters
        ----------
        
        altitude : float, m
        
        
        Examples
        --------
        
        >>> slf.run_max_SAR(altitude)
        performance data
        =========
        LD              = 16.13211    
        turnRadius      = 0.00000     
        TSFC            = 0.00041     
        density         = 0.41270     
        fuelFlow        = 0.19257     
        loadFactor      = 1.00000     
        drag            = 4654.67845  
        lift            = 75089.77219 
        SAR             = 1240.79966  
        thrust          = 4654.67845  
        velocity        = 238.93994   
        EAS             = 0.00000     
        turnRate        = 0.00000     
        climbRate       = 0.00000     
        climbAngle      = 0.00000     
        bankAngle       = 0.00000     
        Mach            = 0.79791     
        """
        atm = ISAtmosphere(altitude)
        Vstall = self.get_Vstall()
        args = (altitude,atm.soundSpeed,)
        bound = np.array([Vstall,0.85*atm.soundSpeed])
        opts = {'maxiter':100,'disp':False}
        def velocity_eqn(V,altitude,soundSpeed):
            D = self.get_required_thrust(V,altitude)
            sfc = self.tm.get_sfc(V/soundSpeed,altitude,D)
            ff = D*sfc #fuel flow kg/sec
            return -V/ff
        rslt = minimize_scalar(velocity_eqn,bracket=bound,bounds=bound,method='Bounded',options=opts,args=args)
        V = rslt.x
        self.set_results(V,altitude)
        return self.results
    
    def run_min_fuel(self,altitude):
        """ Maximum endurance condition
        
        Parameters
        ----------
        
        altitude : float,m
        
        Examples
        --------
        
        >>> slf.run_min_fuel(altitude)
        performance data
        =========
        LD              = 15.99632    
        turnRadius      = 0.00000     
        TSFC            = 0.00036     
        density         = 0.41270     
        fuelFlow        = 0.17296     
        loadFactor      = 1.00000     
        drag            = 4694.19152  
        lift            = 75089.77219 
        SAR             = 1110.09921  
        thrust          = 4694.19152  
        velocity        = 192.00386   
        EAS             = 0.00000     
        turnRate        = 0.00000     
        climbRate       = 0.00000     
        climbAngle      = 0.00000     
        bankAngle       = 0.00000     
        Mach            = 0.64117     
        """
        atm = ISAtmosphere(altitude)
        Vstall = self.get_Vstall()
        args = (altitude,atm.soundSpeed,)
        bound = np.array([Vstall,0.85*atm.soundSpeed])
        opts = {'maxiter':100,'disp':False}
        def velocity_eqn(V,altitude,soundSpeed):
            D = self.get_required_thrust(V,altitude)
            sfc = self.tm.get_sfc(V/soundSpeed,altitude,D)
            ff = sfc*D
            return ff
        rslt = minimize_scalar(velocity_eqn,bracket=bound,bounds=bound,method='Bounded',options=opts,args=args)
        V = rslt.x
        self.set_results(V,altitude) #FIXME: results are not satisfactory
        return self.results
    
    def set_results(self,velocity,altitude):
        results = PerformanceResults()
        fc = FlightConditions(velocity,altitude)
        results.velocity           = fc.velocity
        results.Mach               = fc.Mach
        results.altitude           = altitude
        results.density            = fc.atm.density
        results.thrust             = self.get_required_thrust(velocity,altitude)
        results.drag               = results.thrust
        results.TSFC               = self.tm.get_sfc(results.Mach,altitude,results.thrust)
        results.fuelFlow           = results.TSFC*results.thrust/constants.GRAVITY_ACCEL
        results.SAR                = velocity/results.fuelFlow
        results.lift               = self.bi.wt
        results.LD                 = results.lift/results.drag
        results.weight             = self.ac.get_mass()
        results.emptyWeight        = self.ac.get_mass_empty()
        results.fuelWeight         = self.ac.get_mass_fuel()        
        results.fuelWeightFraction = results.fuelWeight/results.weight
        results.payloadWeight      = self.ac.mass.payload()
        results.reynoldsNumber     = fc.Re
        results.liftCoefficient    = results.lift/self.ac.wing.area/(0.5*fc.atm.density*fc.velocity*fc.velocity)
        ##results.parasiteDrag = self.ac.get_parasite_drag_fw()   
        ##results.induceedDrag 
            
        self.results = results


class ClimbDescent(FlightMechanics):
    """
    Set of methods for climb and descent flight
    
    Examples
    --------
    
    Load aircraft
    
    >>> ac = aircraft_FW.load('X47B')
    
    Initialize steady level flight class

    >>> CLmax3D = 1.6
    >>> clm = SteadyLevelFlight(ac, CLmax3D)
    """
    def run_max_climb_rate(self,altitude):
        """
        calculates maximum climb rate at given altitude
        
        Parameters
        ----------
        
        altitude : float,m
        
        Examples
        --------
        
        >>> clm.run_max_climb_rate(0)
        performance data
        =========
        LD              = 3.18489     
        turnRadius      = 0.00000     
        TSFC            = 0.00027     
        density         = 1.22505     
        fuelFlow        = 21.53547    
        loadFactor      = 1.00000     
        drag            = 42200.72423 
        lift            = 134404.87293
        SAR             = 0.00000     
        thrust          = 78963.87855 
        velocity        = 172.12809   
        EAS             = 0.00000     
        turnRate        = 0.00000     
        climbRate       = 40.24762    
        climbAngle      = 0.23601     
        bankAngle       = 0.00000     
        Mach            = 0.50583    
        """
        atm = ISAtmosphere(altitude)
        Vstall = self.get_Vstall()
        bound = np.array([2.0*Vstall,0.95*atm.soundSpeed])
        opts = {'maxiter':100,'disp':False}
        args = (altitude,)
        fun = lambda velocity,altitude: -self.run_climb_rate(velocity,altitude).climbRate
        rslt = minimize_scalar(fun,bracket=bound,bounds=bound,method='Bounded',
                               options=opts,args=args)
        velocity = rslt.x
        return self.run_climb_rate(velocity,altitude)
    
    def run_most_economical_climb(self,altitude):
        atm = ISAtmosphere(altitude)
        Vstall = self.get_Vstall()
        bound = np.array([Vstall,0.95*atm.soundSpeed])
        opts = {'maxiter':100,'disp':False}
        args = (altitude,)
        def fun(velocity,altitude):
            rslt = self.run_climb_rate(velocity,altitude)
            return -rslt.climbRate/rslt.fuelFlow
        rslt = minimize_scalar(fun,bracket=bound,bounds=bound,method='Bounded',
                               options=opts,args=args)
        velocity = rslt.x
        rslt = self.run_climb_rate(velocity,altitude)
        return rslt
    
    def get_ceiling(self,RCmin=0.5):
        """ 
        calculate altitude at which rate of climb is equal to required value 
        
        Parameters
        ----------
        
        RCmin : float, m/s
            Required rate of climb.
        """
        RCmin = float(RCmin)
        if RCmin<0.0:
            raise ValueError('Rate of climb should be greater than zero')
        fun = lambda h: self.run_max_climb_rate(h).climbRate-RCmin
        xopt = fsolve(fun,5.0e3,xtol=1e-3)[0]
        return xopt

    def get_service_ceiling(self):
        """
        calculate service ceiling
        
        Examples
        --------
        
        >>> clm.get_service_ceiling()
        10084.659906018382
        """
        return self.get_ceiling(0.5)
    
    def get_absolute_ceiling(self):
        """
        calculate absolute ceiling
        
        Examples
        --------
        
        >>> clm.get_absolute_ceiling()
        10221.724848997959
        """
        return self.get_ceiling(0.0)

    def run_max_climb_angle(self):
        pass

    def run_min_glide_angle(self):
        pass
    def run_min_glide_sinkrate(self):
        pass
    def run_climb_rate(self,airspeed,altitude,thrustSetting=1.0):
        basicInput   = self.bi
        thrustModule = self.tm
        fc = FlightConditions(airspeed,altitude)
        V       =fc.velocity
        Cd0     =basicInput.Cd0
        g       =basicInput.g
        k       =basicInput.k
        S       =basicInput.refArea
        W       =basicInput.mass*g
        Q       =fc.dynamicPressure
        tol     =1.0e-4
        change  =tol+1.0
        CA      =0.0 #climb angle, rad
        iMax    =100
        i       =0
        T = thrustSetting*thrustModule.get_thrust_available(altitude)
        while change>tol and i<=iMax:
            L  =W*np.cos(CA)-T*np.sin(CA)
            Cl =L/Q/S
            Cd =Cd0+k*Cl*Cl
            D  =Q*S*Cd
            RC =(T-D)*V/W
            ca =CA
#            if RC/V>1.0:
#                print 'value of RC is too high, temporarily changing the value'
#                print 'this feature MUST be fixed'
#                RC=V#FIXME: must be repaired
            CA =np.arcsin(RC/V)
            change=np.abs(ca-CA)
            i+=1
        results = PerformanceResults()
        results.altitude = altitude
        results.velocity = fc.velocity
        results.thrust = T
        results.climbAngle = np.degrees(CA)
        results.climbRate = RC
        results.Mach = fc.Mach
        results.density = fc.atm.density
        results.lift = L
        results.drag = D
        results.LD = L/D
        if hasattr(altitude,'__iter__'):
            altitude = altitude[0]
        results.TSFC = self.tm.get_sfc(results.Mach, altitude, results.thrust)
        results.fuelFlow = self.tm.get_fuel_flow(results.Mach, altitude, results.thrust)
        results.SAR = results.fuelFlow / results.velocity
        return results


class TurningFlight:
    def run_turn_rate(self,phi,airspeed,loadFactorMax,powerSetting):
        pass


class Field:
    pass


def get_header(string,underline='-'):
    return str(string) + underline*len(string) + '\n'

def get_maxTAS_basic_report(data):
    out   = '{:<23} = {:<12.1f}{:<6}\n'.format('altitude',data.altitude,'m')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('total mass',data.weight,'kg')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('fuel mass fraction',data.fuelWeightFraction*100,'%')
    out  += '{:<23} = {:<12.2f}{:<6}\n'.format('velocity',data.velocity,'m/sec')
    out  += '{:<23} = {:<12.4f}{:<6}\n'.format('Mach',data.Mach,'')
    return out

def get_maxSAR_basic_report(data):
    out   = '{:<23} = {:<12.1f}{:<6}\n'.format('altitude',data.altitude,'m')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('total mass',data.weight,'kg')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('fuel mass fraction',data.fuelWeightFraction*100,'%')
    out  += '{:<23} = {:<12.2f}{:<6}\n'.format('velocity',data.velocity,'m/sec')
    out  += '{:<23} = {:<12.4f}{:<6}\n'.format('Mach',data.Mach,'')
    out  += '{:<23} = {:<12.2f}{:<6}\n'.format('Specific range',data.SAR,'m/kg')
    return out

def get_min_fuel_basic_report(data):
    out   = '{:<23} = {:<12.1f}{:<6}\n'.format('altitude',data.altitude,'m')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('total mass',data.weight,'kg')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('fuel mass fraction',data.fuelWeightFraction*100,'%')
    out  += '{:<23} = {:<12.2f}{:<6}\n'.format('velocity',data.velocity,'m/sec')
    out  += '{:<23} = {:<12.4f}{:<6}\n'.format('Mach',data.Mach,'')
    out  += '{:<23} = {:<12.4f}{:<6}\n'.format('fuel flow',data.fuelFlow,'kg/s')
    return out

def get_max_climb_rate_basic_report(data):
    out   = '{:<23} = {:<12.1f}{:<6}\n'.format('altitude',data.altitude,'m')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('total mass',data.weight,'kg')
    out  += '{:<23} = {:<12.1f}{:<6}\n'.format('fuel mass fraction',data.fuelWeightFraction*100,'%')
    out  += '{:<23} = {:<12.4f}{:<6}\n'.format('rate of climb',data.climbRate,'m/s')
    return out

def point_performance_analysis(ac, altitude=None, filePath=None,
                               fuelFraction=0.5,climbFuelFraction=1.0,detailedReport=False):

    import matplotlib.pyplot as plt
    
    if altitude==None:
        altitude = ac.designGoals.cruiseAltitude
    altitudeRC = 0.0

    CLmaxClean = ac.get_CLmax3D(0)
    ac.get_mission_fuel()
    ac.set_fuel_fraction(fuelFraction)
    

    ## SteadyLevelFlight
    slf = SteadyLevelFlight(ac)
    out = 'POINT PERFORMANCE\n'+'='*40 + '\n'
    out += 'MAXIMUM SPEED'
    out += '\n' + 40*'-' + '\n'
    maxTAS = slf.run_max_TAS(altitude)
    if detailedReport:
        out += maxTAS.__repr__()
    else:
        out += get_maxTAS_basic_report(maxTAS)
    out += '\n'
    out += 'MAXIMUM SPECIFIC RANGE'
    out += '\n' + 40*'-' + '\n'
    maxSAR = slf.run_max_SAR(altitude)
    if detailedReport:
        out += maxSAR.__repr__()
    else:
        out += get_maxSAR_basic_report(maxSAR)
    out += '\n'
    out += 'MINIMUM FUEL FLOW'
    out += '\n' + 40*'-' + '\n'
    minFuel = slf.run_min_fuel(altitude)
    if detailedReport:
        out += minFuel.__repr__()
    else:
        out += get_min_fuel_basic_report(minFuel)
    out += '\n'

    ## ClimbDescent
    ac.set_fuel_fraction(climbFuelFraction)
    clm = ClimbDescent(ac)
    out += 'MAXIMUM CLIMB RATE'
    out += '\n' + 40*'-' + '\n'
    maxRC = clm.run_max_climb_rate(altitudeRC)
    if detailedReport:
        out += maxRC.__repr__()
    else:
        out += get_max_climb_rate_basic_report(maxRC)
    out += '\n'
    serviceCeiling = clm.get_service_ceiling()
    absoluteCeiling = clm.get_absolute_ceiling()
    acRep = clm.run_max_climb_rate(absoluteCeiling)
    
    out += '{:<23} = {:<12.1f}{:<6}\n'.format('Service ceiling',serviceCeiling,'m')
    out += '{:<23} = {:<12.1f}{:<6}\n'.format('Absolute ceiling',absoluteCeiling,'m')
    
    # === operating envelope ===
    nSeg = 20
    S = ac.wing.area
    W = ac.get_mass()
    alt = np.linspace(0,15000,nSeg)
    
    # --- stall speed limit ---
    Mach = np.zeros(len(alt))
    rho = np.zeros(len(alt))
    for i,h in enumerate(alt):
        atm = ISAtmosphere(h)
        rho = atm.density
        Vstall = (2.*W/(rho*S*CLmaxClean))**0.5
        Mach[i] = Vstall/atm.soundSpeed
    # --- maximum speed limit ---
    alt1 = np.linspace(0,acRep.altitude,nSeg)
    Mach1 = np.zeros(nSeg)
    for i,h in enumerate(alt1):
        Mach1[i] = slf.run_max_TAS(h).Mach

    
    plt.figure()
    plt.grid(True)
    plt.hold(True)
    plt.plot(Mach,alt,'ro-')
    plt.plot(acRep.Mach,acRep.altitude,'ks')
    plt.plot(Mach1,alt1,'gs--')
    plt.axis([0,1,0,3e4])
    plt.show()
    # ---
    if filePath==None:
        print out
    else:
        fid = open(filePath,'wt')
        fid.write(out)
        fid.close()
    return out

def run_test1():
    import aircraft_FW
    ac = aircraft_FW.load('X47B')
    print ac.designGoals.fuelMass
    print ac.get_mass_fuel()
    ac.get_mission_fuel()
    print ac.designGoals.fuelMass
    print ac.get_mass_fuel()
    #ac.display()
    slf =SteadyLevelFlight(ac)
    #alt = ac.designGoals.cruiseAltitude
    alt = 1e4

    print 'MAXIMUM SPEED'
    print slf.run_max_TAS(alt)
    print 'MINIMUM SPEED'
    print slf.run_min_TAS(alt)
    print 'MAXIMUM SAR'
    print slf.run_max_SAR(alt)
    print 'MINIMUM FUEL'
    print slf.run_min_fuel(alt)
    
    ac.mass.fuel.set_fuel_fraction(0.5)
    clm = ClimbDescent(ac)
    print 'MAXIMUM CLIMB RATE'
    print clm.run_max_climb_rate(0)
    print 'MOST ECONOMICAL CLIMB'
    print clm.run_most_economical_climb(0)
    print 'SERVICE CEILING'
    print clm.get_service_ceiling()
    print 'ABSOLUTE CEILING'
    print clm.get_absolute_ceiling()


if __name__=="__main__":
    import aircraft_FW
    ac = aircraft_FW.load('X47B')
    point_performance_analysis(ac,0.0, None)